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Abstract 

The binding energies of two-dimensional clusters (puddles) of 4 He are calculated in the framework 
of the diffusion Monte Carlo method. The results are well fitted by a mass formula in powers of 
x = iV -1 / 2 , where N is the number of particles. The analysis of the mass formula allows for the 
extraction of the line tension, which turns out to be 0.121 K/A. Sizes and density profiles of the 
puddles are also reported. 
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I. INTRODUCTION 



In recent years, a great deal of work has been devoted to study quantum liquids in 
restricted geometries.— One important feature of these systems is that their internal structure 
becomes more easily observable than in bulk liquids due to the restricted motion of the 
particles in the confining potential. Among these systems the study of quantum films has 
received particular attention. They consist of liquid helium adsorbed to a more-or-less 
attractive flat surface. In 1973 M. Bretz et al£ observed for the first time the adsorption 
of 4 He onto the basal plane of graphite. In the last few years, adsorption properties of 
helium on different substrates such as carbon, alkali and alkaline-earth flat surfaces, carbon 
nanotubes and aerogels have become a fertile topic of research. 

The structure and growth of thin films of 4 He adsorbed to a substrate was studied by 
Clements et a/.— employing the optimized hypernetted-chain Euler-Lagrange theory with 
realistic atom-atom interactions. It turns out that films with low surface coverages, where all 
atoms cover the surface with a thickness corresponding to a single atom, can be approximated 
reasonably well by a 2D model. In connection with these systems, an interesting question 
naturally arises as how physics depends on the dimensionality of the space. 

The homogeneous 2D liquid has been studied using different theoretical methods, such 
as molecular dynamics 4 - and quantum Monte Carlo simulations either Green's Function^ 
or diffusion^ techniques. The inhomogeneous case was studied by Krishnamachari and 
Chester,^ who used a shadow variational wave function to describe 2D puddles of liquid 
4 He. In this work we report energies and density profiles of puddles calculated within the 
diffusion Monte Carlo (DMC) method. Our main objective is to give an accurate estimate 
of the line energy or the line tension of the 2D liquid 4 He. As atom-atom interaction we 
have used the revised version of the Aziz potential dubbed as HFD-B(HE).— This potential 
has been used to study ground-state properties of 3D bulk 4 He and 3 He,^*i£ within the DMC 
framework, and it has proved to accurately reproduce the ground-state properties of both 
liquids at zero temperature. 

The trial wave function used for the importance sampling in the DMC calculation is 
introduced in Section II, where the variational Monte Carlo results for this wave function 
are also reported. A brief explanation of the DMC techniques used in the present paper is 
presented in Section III. Section IV contains the DMC results and their analysis in terms 
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of a mass formula in 2D. The line tension is extracted from this mass formula. Properties 
characterizing the puddles, such as the density profiles, are discussed in Section V. Finally, 
the main conclusions are summarized in Section VI. 



II. VARIATIONAL GROUND-STATE ENERGIES 

To study a system of N 4 He atoms in two dimensions we start from the following trial 
wave function 



written in the same way as in the 3D case.— The coordinate R indicates the set of coordinates 
of all the particles {ri, r 2 , r^}, while stands for the interparticle distance, = |rj— r,[ . 
The trial wave function contains the simple McMillan form^ 2 - to deal with the very short- 
range part of the interaction, and the translationally invariant part of a harmonic oscillator 
(HO) wave function with parameter a, to roughly confine the system. 

In our calculations the value h 2 /m 4 = 12.1194 K A 2 has been employed for the atom 
mass and the parameters b and v have been fixed to the values 3.00 A and 5, respectively, 
the same values as in 3D calculations. The variational search has thus been restricted to 
the HO parameter a, whose optimal value is given in Table ID The expectation value of 
the Hamiltonian, as well as the separate contributions of kinetic and potential energies are 
given in the same Table for puddles with N atoms. It can be seen that the total energy 
results from an important cancellation between kinetic and potential energies, which is in 
fact larger than in the 3D case. Let us recall that in 3D bulk, the energy per particle results 
from adding « 14 K of kinetic energy with « —21 K of potential energy. In 2D, both kinetic 
and potential contributions are very close to each other, which makes the calculation very 
delicate. 

In the last column of Table |U the VMC results of Krishnamachari and Chested are 
reported. As compared with their results, our calculations provide smaller binding energies, 
in spite of the fact that the interaction used in Ref. His an older version of the Aziz potential, 
which tends to underbind the systems. In fact, the shadow wave function used in Ref. [3 
contains more elaborated correlations not present in our simple trial wave function. The 
VMC energy for the bulk system corresponds to the saturation density, po = 0.04344 A -2 , 
taken from the DMC calculation of Ref. la. 




(1) 
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We have also performed calculations using a different trial wave function, replacing the 
translationally invariant gaussian part by an exponential one, i.e. 

I ( b_Y a 



$ T (R) = IJexp 

i<3 



2 In J 



(2) 



expecting that this larger tail in the wave function will result in more binding. Actually, we 
do not find significant differences for small values of N. For instance, in the case N = 8, 
using the same values for b and v as before, we get E/N = —0.2178(5) K, T/N = 1.266(2) 
K, V/N = -1.484(2) K for a = 0.035 A" 1 . When the values of 6, v and a are optimized, 
we obtain a slightly larger binding energy, E/N = —0.2267(8) K for b = 3.04 A, u — 5.0 
and a = 0.035 A -1 . For greater values of N, the gaussian ansatz tends to provide more 
binding than the exponential. For example, with the exponential ansatz, for iV = 16 we get 
E/N = —0.1816(7) K for a = 0.023 A -1 , and optimizing the different parameters one gets 
E/N = -0.2514(6) K, with b = 3.04 A. In conclusion, the gaussian wave function seems 
appropriate to be used as importance sampling in the DMC calculations. 



III. DIFFUSION MONTE CARLO GROUND STATE ENERGIES 

Quantum Monte Carlo (QMC) methods provide the exact ground-state energy of a boson 
system, except for statistical errors. These techniques solve numerically the Schrodinger 
equation by means of a statistical simulation. They have been widely described in the 
literature; hence we briefly recall here the main ideas, referring the reader to, for example, 
Ref.Qfor a more detailed description on QMC techniques. In this work we use the diffusion 
Monte Carlo (DMC) method to solve the Schrodinger equation in imaginary time (r = it) 
for the function 

/(R,r) = $ r (R)*(R,T) (3) 

where R represents all the particle coordinates and is usually called "the walker", \I/(R, r) is 
the wave function of the system, and $^(R) is the previously determined trial wave function 
(Sect. II), used here as importance sampling. It is convenient to write the solution of the 
time-dependent Schrodinger equation in the following form 

/(R, r + At) = J dR'G(R, R', Ar)/(R', r) (4) 
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where G is the time-dependent Green's function and is formally written as 

G(R, R'; r) = (R | e~ Hr | R') (5) 

where H is the Hamiltonian of the system. The function G(R, R'; r) represents the ampli- 
tude probability for the transition from an initial state R' to a final one R in a time r. In 
the limit r — > oo, Eq. (jlj) gives the exact ground-state wave function. Thus, knowing G for 
infinitesimal time steps At, the asymptotic solution for large times /(R, r — > oo) can be 
obtained by solving iteratively the above equation. To this end, the exponential entering 
Eq. © is approximated to some fixed order in Ar. Both first and second ordei^ propagators 
have been implemented in the present work and both of them provide the same extrapolated 
energy, within statistical errors, using the trial function $t of Eq. (jlj) as guiding function. 
Our simulations have been carried out with a population of typically 400 walkers. As usual, 
some runs are first done to establish the asymptotic region of the short time propagator, then 
several values of the time step have been used, and finally a fit, either linear or quadratic, 
has been carried out to obtain the extrapolated energy. For example, for N = 16 the time 
steps 0.0001, 0.0002, 0.0003, and 0.0004 have been used to perform the extrapolation. In 
general, the statistical error is of the order of the time step error in our calculations. 

In Table HU we present the results of our linear DMC calculations of the total energy per 
particle for puddles containing N atoms. We have also reported and reproduced the results 
of the binding energy per particle of homogeneous 2D liquid 4 He at the equilibrium density 
pDMC _ 0.04344(2) A -2 ; obtained in Ref. |6j, where the same version of the Aziz potential 
was used. For this case, the simulations have been carried out for a system of 64 atoms 
with periodic boundary conditions, for which the errors due to finite size effects are smaller 
than the statistical ones^. We have also performed quadratic DMC calculations for some 
puddles, and found results which are compatible with the linear DMC ones within their 
error bars. For example, the quadratic algorithm provides E quad = —0.2612(2) K for = 8 
and —0.652(4) for N = 64. As expected, the DMC results lower the corresponding energies 
obtained by VMC either with our simple variational wave function or with a shadow wave 
function,— by up to ~ 25% in the case of the bulk system. It is worth mentioning that the 
final DMC result for the energy does not depend on the trial wave function for a boson 
system like the studied here, a fact that in the present case has been numerically checked 
for the Gaussian and the exponential ansatzs, Eqs. (|XJ) and (j2j). Indeed, for boson systems 
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the DMC method provides exact ground-state energies, within statistical errors. 

IV. ENERGY AND LINE TENSION 

For a saturating self-bound system, the ground-state energy per particle can be expanded 
in a series of powers of the variable iV _1//D , where N is the number of constituents and D is 
the dimensionality of the space. This is the well-known mass formula, which in the present 
case writes 

E(N)/N = e b + e l x + e c x 2 + ••• (6) 

with x = iV~ 1//2 . The two first coefficients of this expansion are the bulk energy e b and the 
line energy £j, out of which the line tension A is defined by 27rr A = £j. Here r is the unit 
radius, defined as the radius of a disk whose surface is equal to the inverse of the equilibrium 
density of the 2D bulk liquid, i.e. po^r^ = 1. Finally, e c is the so-called curvature energy. 

Our calculated ground-state energies (Tables |U and HH) are plotted in Fig. [T] as a function 
ofiV- 1 / 2 . One can see that the differences between our VMC and DMC energies increase 
with the number of atoms in the puddle. This clearly mirrors the simplicity of the trial 
wave function, which could be improved by including, for example, three-body correlations. 
Nevertheless this trial function is adequate for the importance sampling in the DMC calcu- 
lation. 

We have fitted these energies to a parabolic mass formula like Eq. JBJ). The coefficients 
of the fit are given in Table together with the deduced line tension. Notice that the 
coefficient e b is identical, within statistical errors, to the bulk energy per particle of Table ITT1 
In fact, the \ 2 °f the fit is very small, \ 2 — 5-7 x 1CT 6 . Regarding the line tension, and 
despite using a different version of the Aziz potential and a different trial function, we notice 
that our VMC estimate is rather close to the one reported in Ref. 7 , A = 0.07 K/A. However, 
both VMC results are remarkably different from the DMC line tension, A = 0.121 K/A. 

To stress the curvature effect we have also plotted in the figure a straight line between 
the N = 8 and bulk DMC values. A linear fit of the DMC energies gives coefficients 
Sb = —0.885 K and E\ = 1.80 K, which are appreciably different from the previous ones. The 
bulk energy extrapolated from this linear fit differs from the directly calculated value, and 
the corresponding line energy is closer to the variational one, thus giving a bad estimation 
for the linear tension. In all cases, the line energy coefficient is approximately minus twice 
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the volume energy term, similarly to the 3D case, 11 and therefore one expects curvature 
effects to be important. In both VMC and DMC cases the extracted e c is negative, i.e. the 
binding energy is a convex function of x as it also happens for the 3D clusters^. This is 
in contrast with the value of e c reported in Ref. Q which was positive but rather smaller in 
absolute value and with larger error bars. Actually, as argued in Ref. [ll| for 3D clusters, 
one would expect the curvature correction to the energy of a circular 2D cluster to be 
positive. Therefore, one should take the extracted value for e c with certain caution and not 
to emphasize its physical significance. However, it turns out that the value and sign of e c 
are stable against different possible fits, e.g., changing the number of points to build the fit, 
or using a cubic mass formula. In any case the two first coefficients £&, E\ are quite robust 
against all performed fits. As an illustration, if one takes out the bulk point, the predicted 
bulk energy per particle and surface tension are equal to the reported ones within 1% and 
5% respectively. Therefore, the extracted line tension should be reliable, as it also happens 
for 3D clusters^. 



V. DENSITY PROFILES 



The calculation of observables given by operators that do not commute with the Hamilto- 
nian poses new problem to the DMC method. After convergence, the walkers are distributed 
according to the so-called mixed probability distribution given by the product of the exact 
and the trial wave functions. Therefore averaging the local values of the operator does not 
give the exact expectation value unless the operator commutes with the propagator. The 
result obtained by straightforward averaging is the mixed estimator which is of first order 
error in the trial wave function. Several options have been proposed in the literature in 
order to obtain unbiased (trial function independent and exact) values. In this work we 
have used the so called forward or future walking technique^ to calculate unbiased, also 
called pure, density profiles. The key ingredient to correct the mixed estimator is to include 
as a weight in the sampling the quotient $ exact (R)/$triai(R) for each walker, given by the 
asymptotic number of walkers. Several algorithms have been proposed in order to compute 
this quantity. In this work we use the algorithm developed in Ref. Q that constitutes a 
simple and efficient implementation of the future walking method. 

The pure DMC estimates of the density profiles for several puddles are plotted in Fig. [21 
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The figure also contains an horizontal line which indicates the saturation density (Po MC ) of 
the homogeneuos system. For the puddle containing 36 atoms, the VMC profile obtained 
from a Gaussian ansatz (Eq. (JTJ) ) is also shown for comparison as a dotted line. As one can 
appreciate in the figure, the process of optimization implied by the DMC method changes 
the profile reducing its thickness, i.e., producing a sharper surface. It can be seen that 
for the smaller clusters the central density is below po, while for the larger values of N 
shown in the figure the central density is above po, indicating a leptodermous behaviour. 
One expects that, increasing the number of particles, the central density will approach po 
from above, as in the 3D case . 1 ^ 1 ^ It is worth noticing the oscillating behaviour in the 
interior part of the density profile for N = 121. It is difficult however to decide whether 
these oscillations are genuine or are simply due to a poor statistics in evaluating the pure 
estimator. Unfortunately, to discard this last option would require an exceedingly long 
computing time, within the scheme of this work. 

The solid lines plotted in Fig. [2] are fits to our DMC densities provided by a generalized 
Fermi profile of the form: 

P(r) = T . PJ 7T^ (7) 



1 + exp 

The parameters defining the Fermi profile are given in Table ITVl together with the thickness 
t and the root mean square (rms) radius (r 2 ) 1 / 2 . We have checked that the rms radius 
calculated within the DMC code and the one derived form the fit agree to better than 0.5%, 
except for the N = 121 case, where the difference is 1%. The rms radius grows with the 
number of particles as iV 1 / 2 , as expected. Therefore it grows faster than in 3D, in which 
case grows as iV 1 / 3 . This behaviour allows for an alternative determination of the saturation 
density by performing a linear fit to the relation 



N 

The extracted value of p from the mean square radius reported in Table Hvl is 0.043 A -2 , in 
good agreement with the determination from the calculation for the homogeneous system. 

In the interval of N considered, the thickness t, defined as the distance over which the 
density falls from 0.9 of its value at origin to 0.1, is continuously increasing. However, as the 
finite value of the thickness for the semiinfinite system should define the asymptotic value 
of t, one expects that for larger puddles the thickness will probably have a maximum and 
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smoothly approach this asymptotic value, as happens in the 3D case.— Finally, one also 
observes the asymmetric character of the density profiles with respect to the point at which 
the density falls at half its value at the origin. This can be appreciated by looking at the 
value of v, which grows with N, and also to the increasing difference between the quantities 
R and (r 2 ) 1 / 2 . 

VI. SUMMARY AND CONCLUSIONS 

In this work we have considered strictly two-dimensional systems of liquid 4 He, which are 
of course an idealization of a real quantum film. They are nevertheless interesting because 
their study can enlighten the underlying structure of real quasi-2D systems. Of course, in 
the latter case, one has to take also into account the interaction with the substrate, which 
basically provides a global attractive potential. In the ideal 2D case, the suppression of 
the wave function component in the third dimension, produces an increment of the global 
repulsion between atoms, resulting in a smaller binding energy per particle, and a decrease 
of the equilibrium density.— 

The binding energies of two-dimensional 4 He clusters, calculated by means of a diffusion 
Monte Carlo method, are well fitted by a mass formula in powers of x = N~ l l 2 . The 
analysis of the mass formula provides the main result of this paper, namely the value of the 
line tension A = 0.121 K/A, which significantly differs from the one obtained from a similar 
analysis of VMC data and the one previously reported in the literature^. The quadratic term 
of the mass formula cannot be neglected and results in a negative value of the curvature 
energy as in the 3D case.— J£ However, the studied clusters may be too small to give physical 
significance to this result. 

The density profiles obtained with the pure estimator have been fitted by a generalized 
Fermi function, and the behaviour of the rms radius and the thickness as well as the asym- 
metry character of the profile as a function of N have been discussed. However, calculations 
for larger puddles, which are out of the scope of the present paper due to limitations in 
computing time, would be necessary to describe the complete iV-dependence of the density 
profiles. 
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TABLE I: Variational results for the ground-state energy per particle E/N of 2D 4 He puddles of 
various cluster sizes. The confining HO parameter a is given in A -1 and all energies are in K. The 
expectation values of the kinetic and the potential energies are also displayed. The column labelled 
KC refers to the VMC results of Ref . U . 



N 


a 


E/N 


T/N 


V/N 


KC 


8 


0.1565 


-0.2239(2) 


1.3003(6) 


-1.5242(5) 




16 


0.129 


-0.3510(2) 


1.7354(6) 


-2.0864(5) 


-0.380(8) 


36 


0.094 


-0.4532(4) 


2.031(3) 


-2.484(3) 


-0.471(7) 


64 


0.073 


-0.4961(7) 


2.159(2) 


-2.655(2) 


-0.528(5) 


121 


0.054 


-0.5241(6) 


2.223(2) 


-2.747(2) 


-0.570(7) 


165 


0.047 


-0.5328(3) 


2.289(1) 


-2.822(1) 


-0.602(7) 


512 0.0266 


-0.5493(5) 


2.282(3) 


-2.831(3) 


-0.621(2) 


oo 


0.0000 


-0.6904(8) 


4.312(2) 


-5.003(1) 





TABLE II: Energy per particle (in K) for 2D 4 He puddles for various cluster sizes obtained with 
the DMC algorithm. 



N 
E/N 



8 16 36 64 121 oo 

-0.2613(4) -0.4263(4) -0.578(2) -0.658(4) -0.710(2) -0.899(2) 



TABLE III: Coefficients (in K) of a parabolic fit of the mass formula, as given in Eq. The last 
column displays the deduced line tension (in K A -1 ). 





£6 £l 




A 


VMC 


-0.654(1) 1.41(1) 


-0.62(2) 


0.083(1) 


DMC 


-0.898(2) 2.05(2) 


-0.71(3) 


0.121(1) 
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TABLE IV: Parameters of a Fermi-profile fit to the density profiles. All lengths are in A and pf is 
in A -2 . The parameter v is adimensional. 



N 


pf R c v 


t 




8 


0.03740 9.308 2.156 1.739 


8.166 


7.20 


16 


0.04204 13.38 2.656 2.284 


9.580 


9.18 


36 


0.04305 19.47 3.104 2.400 


11.11 


12.91 


64 


0.04386 26.68 3.783 3.111 


13.09 


16.68 


121 


0.04304 40.09 5.566 4.714 


18.52 


23.15 




J I I I I I I I I 

0.1 0.2 0.3 0.4 



FIG. 1: Energies per particle (in K) of iV-atom puddles as a function of N~ 1 / 2 , obtained from our 
VMC (squares) and DMC (circles) calculations. The interaction used is Aziz HFD-B(HE). Dashed 
and solid lines correspond to a least square fit to these energies. The dot-dashed line is a straight 
line between the N = 8 and bulk DMC values. 
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r(A) 

FIG. 2: Density profiles for 4 He puddles with various number of atoms, N = 8 (circles), 16 (squares), 
36 (diamonds), 64 (triangles up) and 121 (triangles down), obtained from our pure estimators 
for the linear DMC calculations. The solid horizontal line indicates the saturation density of the 
homogeneous system. The dotted line is the VMC profile for N = 36 with a Gaussian trial function. 
The figure also contains the fits to the data provided by a generalized Fermi function, as explained 
in the text. 
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